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HIGHLIGHTS 


•  Improved  method  to  determine  averaged  transport  parameters  e.g.  porosity,  tortuosity. 

•  Based  on  reconstruction  of  3D  micro-structures  of  a  metal-supported  SOFC. 

•  Correct  lateral  boundary  conditions  change  state-of-the  art  approach  with  47%. 

•  Assessment  of  anisotropic  micro-structures  by  determining  the  porosity  distribution. 

•  Transport  parameters  for  different  dense  layers  are  computed  separately. 
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Metallic  supported  Solid  Oxide  Fuel  Cells  (SOFCs)  are  considered  as  a  durable  and  cost  effective  alter¬ 
native  to  the  state-of-the-art  ceramic  supported  cell  designs.  In  order  to  understand  the  mass  and  charge 
transport  in  the  metal-support  of  this  new  type  of  cell  a  novel  technique  involving  X-ray  tomography  and 
micro-structural  modelling  is  presented  in  this  work.  The  simulation  technique  comprises  a  novel 
treatment  of  the  boundary  conditions,  which  leads  to  more  accurate  effective  transport  parameters 
compared  to  those,  which  can  be  achieved  with  the  conventional  homogenisation  procedures. 
Furthermore,  the  porosity  distribution  in  the  metal-support  was  determined,  which  provided  informa¬ 
tion  about  the  inhomogeneous  nature  of  the  material.  In  addition  to  that,  transport  parameters  for  two 
identified,  different  dense  layers  of  the  metal-support  are  evaluated  separately. 

The  results  of  the  evaluation  show  three  significant  findings.  Firstly,  that  the  effective  transport  pa¬ 
rameters  are  up  to  a  factor  of  20  lower  compared  to  other  SOFC  anodes.  Secondly,  that  the  micro¬ 
structure  can  have  regions,  which  are  much  denser,  especially  the  first  100  pm  of  the  interface  be¬ 
tween  gas-channel  and  support-structure.  Thirdly,  that  the  calculation  of  the  transport  parameters  de¬ 
pends  on  the  correct  application  of  boundary  conditions. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  solid  oxide  fuel  cell  (SOFC)  technology  offers  effective 
conversion  of  chemical  stored  energy  to  electricity.  The  state-of- 
the-art  SOFC-stack  is  based  mainly  on  all-ceramic  cells,  which  are 
brittle  and  thus  challenging  to  operate.  During  the  operation  of  the 
SOFCs  heat  is  generated,  which  must  be  transported  out  of  the  stack 
resulting  in  thermal  gradients.  These  thermal  gradients  will  cause 
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thermal  stresses,  which  may  induce  cracks  in  the  ceramic  cells.  A 
second  failure  mechanism  in  state-of-the-art  cells  is  the  re¬ 
oxidation  of  the  nickel/zirconia  cermet  anode  and  substrate.  A  re¬ 
oxidation  of  the  nickel  will  result  in  a  volume  expansion  and  thus 
the  formation  of  cracks.  To  overcome  these  drawbacks,  alternative 
types  of  SOFCs  based  on  porous  metallic  supports  have  been 
developed,  exhibiting  a  number  of  advantages  as  an  improved 
mechanical  strength,  a  higher  thermal  conductivity  of  the  substrate 
and  a  higher  re-oxidation  tolerance  of  the  applied  alloys.  The  actual 
cell  and  stack  developments  performed  worldwide  have  shown 
very  promising  results  compared  to  the  state-of-the-art  cells  [1—9]. 

In  order  to  avoid  failures  and  optimise  the  cells  further,  exper¬ 
iments  are  carried  out  to  study  their  behaviour  [1,10],  but  the 
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theoretical  foundation  to  understand  their  behaviour  must  also  be 
developed.  To  this  end,  multi-physics  models  describing  the  variety 
of  physics  are  being  developed  see  e.g.  Refs.  11  — 13  .  The  physics 
involved  are:  electronic  and  ionic  charge  transfer;  solid  mechanics; 
heat  transfer;  gas  diffusion  and  flow.  The  parameters  used  in  these 
multi-physics  models  can  either  be  measured  or  deducted  by 
micro-structural  modelling.  The  measurements  are  typical  on  the 
cell  level,  why  a  convoluted  response  is  measured  and  the  physi¬ 
cally  meaningful  parameters  can  only  be  determined  by  a  reverse 
engineering  approach  if  an  appropriate  electro-chemical  model  is 
existing  [14  .  In  addition  to  that,  if  a  component  is  modified  or  a 
new  component  is  developed  all  physical  parameters  have  to  be  re¬ 
measured.  Therefore,  a  novel  route  based  on  micro-structural 
modelling  of  the  material  parameters  is  pursued  in  this  work. 

The  various  physical  phenomena  occurring  in  a  SOFC  take  place 
on  a  large  span  of  length  scales  from  sub-micron  to  centimetres. 
Thus,  in  order  to  obtain  meaningful  computational  results,  the 
material  properties  and  the  material  structure  of  the  SOFC  under 
investigation  have  to  be  represented  in  detailed  micro-structural 
models.  Flowever,  a  high  spatial  resolution  of  the  electrode  struc¬ 
ture  will  result  in  a  computational  model  containing  a  tremendous 
amount  of  mesh  cells,  and  will  thus  lead  to  impractical  calculation 
times.  Calculations  of  an  entire  SOFC-stack  (12  cm  x  12  cm  x  10  cm) 
with  a  resolution  of  1  pm3  would  result  in  approx.  1  x  1015 
computational  cells,  which  would  require  memory  capacities 
exceeding  most  of  the  available  supercomputers. 

Often  the  details  occurring  on  smaller  scales  are  indirectly 
included  through  the  applied  material  models  or  so-called  consti¬ 
tutive  models,  in  which  the  material  parameters  are  fitted  to 
experimental  results.  With  the  many  interacting  physics  in  SOFC, 
the  modelling  entails  a  large  number  of  unknown  parameters.  The 
combined  response  of  the  cell  can  to  some  degree  be  decoupled  by 
electrochemical  impedance  spectroscopy  and  e.g.  by  controlling  the 
gas  composition  on  the  anode  and  the  cathode  side,  such  that  pa¬ 
rameters  related  to  phenomena  occurring  either  side  can  be 
determined.  However,  the  combined  response  is  a  convolution  of 
all  the  parameters,  which  makes  the  fitting  task  quite  difficult. 

This  was  successfully  performed  and  shown  in  Ref.  [14]  at 
elevated  temperatures  (>800  °C).  However,  at  lower  temperatures 
the  gas  diffusion  process  in  the  anode  substrate  is  overlapped  by  a 
process  related  to  the  hydrogen  electro-oxidation  in  the  imped¬ 
ance  spectra,  and  deconvolution  is  thus  not  possible.  The  cell 
studied  in  this  work  has  limited  stability  at  these  temperatures, 
and  thus  so  far  it  has  not  been  possible  to  determine  the  mass 
transport  resistance  with  electrochemical  impedance  spectros¬ 
copy.  Alternatively,  one  of  the  processes,  e.g.  the  gas  diffusion, 
could  be  studied  separately  allowing  for  a  better  deconvolution  of 
the  overall  response.  It  should,  however,  be  noted  that  the 
deconvolution  approach  only  allows  for  determining  the  inte¬ 
grated  diffusion  resistance  of  all  layers  on  either  side  of  the  elec¬ 
trolyte.  Thus,  for  a  direct  study  of  the  implication  of  variation  of  the 
composition  and  thickness  of  the  individual  layers  other  methods 
must  be  applied.  Furthermore,  it  is  given  that  the  usage  of  this 
measurement/fitting  approach  will  only  re-produce  the  experi¬ 
ments  that  they  rely  on,  but  the  impact  of  variations  of  micro¬ 
structure  or  maybe  even  layer  thicknesses  requires  new  experi¬ 
ments.  In  the  spirit  of  John  von  Neumann,  a  mathematician  and 
computer  scientist: 

With  four  parameters  I  can  fit  an  elephant,  and  with  five  I  can 

make  him  wiggle  his  trunk; 

computational  attempts  to  reduce  the  level  of  uncertainty  about 
these  parameters  should  be  undertaken. 

The  possibility  to  analyse  the  micro-structure  of  the  elec¬ 
trodes  by  e.g.  focused  ion  beam  -  scanning  electron  microscope 


(FIB-SEM)  or  X-ray  computed  tomography  (CT)  has  therefore  lead 
to  a  number  of  modelling  studies  of  the  phenomena  occurring  on 
the  micro-scale  15-21  ].  The  established  methodology  is  to  apply 
FIB-SEM  measurements  for  regions,  where  a  high  resolution  is 
needed  see  e.g.  Ref.  [22].  Only  a  small  part  of  the  electrode  can  be 
investigated  with  this  method,  whereas  CT-scanning  can  be  used 
to  analyse  larger  volumes  of  the  micro-structure.  Nevertheless, 
X-ray  nano-tomography  yields  comparable  results  to  FIB-SEM 
[23,24].  The  two  approaches  compliment  and  overlap  each 
other  in  their  regime  of  use.  Common  for  both  is  that  the  data 
obtained  is  a  sequence  of  two  dimensional  pictures,  which  has  to 
be  transformed  into  a  three-dimensional  voxelised  picture  of  the 
micro-structure.  This  is  used  to  reconstruct  a  computational  ge¬ 
ometry,  after  image  processing  with  an  in-house  algorithm  based 
on  [25-28  .  Subsequently,  effective  transport  properties  are 
calculated  using  numerical  models.  This  approach  has  been  used 
a  number  of  times  for  SOFC  anodes  and  cathodes  [15-21  ,  but 
was,  to  the  author's  knowledge,  not  yet  applied  to  the  porous 
metal-support.  The  new  technology  of  a  metallic  porous  support 
for  the  SOFC  implicates  a  change  of  the  micro-structure  due  to 
corrosion  during  operation.  This  is  a  very  important  issue  in  the 
development  of  durable  metal-supports  and  is  object  to  many 
experiments  [29-32  .  While  improving  the  stability  of  the  cells,  a 
theoretical  study  of  the  influence  of  the  changing  micro¬ 
structure  on  the  mass  transport  is  needed  to  study  the  feasi¬ 
bility  of  this  new  technology.  This  could  be  done  by  such  a  micro- 
structural  model  which  is  capable  of  determining  effective 
transport  parameter  that  can  be  used  for  multi-physics  models  in 
cell  simulations. 

The  diffusion  through  the  micro-structure  of  the  porous  support 
can  be  determined  by  use  of  a  numerical  method,  e.g.  computa¬ 
tional  fluid  dynamics  (CFD),  Lattice  Boltzmann  method  (LBM)  or 
the  finite  element  method.  However,  it  would  be  unfeasible  with 
current  computational  power  to  represent  a  full  SOFC  stack  for  that 
sake  by  this  approach.  Thus,  in  a  multi-scale  modelling  approach 
the  accurate  micro-structural  models  can  be  utilised  to  determine 
the  response  of  the  smaller  scale  model.  The  results  can  be  repre¬ 
sented  through  homogenised  material  parameters  in  homogenised 
continuum  models  [33].  Hereby  an  accurate  description  of  the  flow 
is  ensured  together  with  a  reduction  of  computational  cells  and 
required  computer  time  by  a  factor  of  105  or  more.  By  this  approach 
more  fundamental  knowledge  about  the  material  parameters  and 
the  occurring  physical  phenomena  are  obtained. 

In  homogenisation  a  boundary  value  problem  is  posed  and 
solved  to  obtain  the  homogenised  material  parameters  [34-36]. 
The  boundary  conditions  for  the  flow  problem  can  either  be  the 
flux,  concentration  or  a  periodic  boundary  condition.  It  has  become 
customary  to  compute  the  effective  diffusion  coefficient  for  diffu¬ 
sion  through  the  thickness  of  the  cell  by  applying  symmetry  con¬ 
dition  (concentration  gradient  perpendicular  to  the  boundary 
equals  zero)  for  the  lateral  boundary  (parallel  to  the  main  diffusion 
path)  [15-19].  This  effectively  introduces  dead-end  pores,  or  in 
other  words  blocks  tortuous  paths  involving  diffusion  in  and  out  of 
the  domain  in  the  vicinity  of  the  lateral  boundary.  The  only  work¬ 
around  so  far  is  to  perform  the  simulation  on  a  sufficiently  large 
volume  that  suppresses  any  impact  of  the  lateral  boundaries,  but 
requires  a  much  larger  volume  to  be  reconstructed  and  more 
computing  power  and  time  to  perform  the  simulation.  In  Ref.  [16  it 
was  shown  that  with  increasing  size  of  the  base  area  the  values  of 
interest  converge  because  the  error  due  to  the  transport  via  the 
lateral  boundaries  decreases.  Nevertheless,  no  comparison  be¬ 
tween  periodic  and  zero  gradient  boundary  conditions  and  their 
influence  was  made.  So  far  the  effect  of  different  boundary  condi¬ 
tions  for  the  top  and  bottom  boundaries  is  analysed  by  Ref.  20],  but 
not  for  the  lateral  faces. 
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In  this  work  homogenised  parameters  of  a  metal-support  of  a 
metallic  supported  SOFC  are  determined.  This  will  be  a  first  step 
toward  a  degradation  model  in  order  to  understand  the  impact  of 
corrosion  on  the  long  term  performance  of  metallic  supported 
SOFCs  by  obtaining  the  effective  diffusion  coefficient  of  the  porous 
structure,  which  is  an  important  input  parameter  for  the  simulation 
of  a  multi-scale  SOFC  model.  The  parameters  are  determined  with  a 
more  accurate  method  for  homogenisation,  where  the  domains  in 
the  vicinity  of  the  lateral  boundaries  are  not  blocked  for  flow,  i.e.  by 
use  of  periodic  boundary  conditions.  Furthermore,  the  influence  of 
the  inaccurate  application  of  boundary  conditions  on  the  effective 
transport  properties  will  be  investigated.  Also  the  concepts  of  tor¬ 
tuosity  r  and  the  tortuosity-factor  r2  are  discussed  in  detail  for 
further  correct  use  of  the  results  determined  in  this  work.  The 
inhomogeneous  nature  of  the  micro-structure  will  be  addressed 
and  analysed  by  determining  the  porosity  distribution.  In  addition 
to  that,  the  transport  parameters  of  the  different  dense  layers  are 
evaluated,  which  is  relevant  in  inhomogeneous  micro-structures. 
Finally,  the  size  of  a  stochastic  equivalent  representative  volume 
element  (SERVE)  [15,16,20]  and  its  limitations  in  inhomogeneous 
micro-structures  are  discussed. 

2.  Reconstruction  of  the  micro-structure 

In  this  paper  the  porous  support  of  a  metal-supported  SOFC 
produced  at  Technical  University  of  Denmark  is  analysed.  The 
support  is  made  by  tape  casting  of  a  slurry  with  metallic  powder 
and  by  subsequent  sintering  in  a  hydrogen  furnace.  For  processing 
details  please  refer  to  Ref.  [37  .  A  comprehensive  description  of  the 
entire  cell  can  be  found  in  Ref.  [2],  a  detailed  characterisation  via 
electro-chemical  impedance  spectroscopy  is  reported  in  Ref.  [38]. 
The  sample  is  a  sintered  half-cell  (porous  metal-support,  anode 
functional  layer  and  electrolyte)  and  has  not  been  exposed  to  any 
operation,  thus  no  corrosion  has  taken  place  in  it. 

The  reconstructed  micro-structure  of  this  metal-support  is  ob¬ 
tained  by  X-ray  CT  scanning  and  by  use  of  an  in-house  algorithm 
[39]  for  generating  a  computational  mesh.  The  imaging  was  done 
with  a  SkyScan  1172  X-ray  micro-CT  system,  equipped  with  a 
100  keV  X-ray  source  and  an  11  megapixel  (4000  x  2670  pixel)  X- 
ray  detector.  A  volume  of  1.1  mm  x  2.2  mm  x  0.37  mm  was  scanned 
with  a  frame  averaging  of  two  and  a  rotation  step  of  0.36,  so  that 
the  entire  thickness  of  the  actual  metal-support  is  covered.  The  X- 
ray  CT  scan  provides  a  stack  of  images  a  single  of  these  can  be  seen 
in  Fig.  1  together  with  its  segmented  counterpart,  where  the  ma¬ 
terial  is  white  and  the  pores  are  black.  The  gas-channel  will  be  on 
top  of  this  sintered  half  cell.  The  lower  part  comprises  the  anode 
functional  layer  (AFL)  and  the  gas  tight  electrolyte.  The  evaluation 


in  this  work  is  limited  to  the  porous  metal-support.  The  in-house 
algorithm  regenerates  the  micro-structure  into  a  computational 
mesh,  by  analysing  the  2-D  grey-scale  images  obtained  with  X-ray 
CT-scans.  These  2-D  images  are  combined  to  a  3-D  grey-scale  im¬ 
age.  The  images  are  digitised  and  reconstructed  using  statistical- 
functions  [26-28].  With  this  information  the  algorithm  de¬ 
termines  positions,  where  pores  or  metal  are  prevailing.  Thus,  a  3-D 
voxel  mesh  can  be  reconstructed,  which  represents  the  real  micro¬ 
structure.  A  3-D  plot  of  the  micro-structure  obtained  by  this 
method  is  shown  in  Fig.  2.  The  obtained  result  is  written  in  a  format 
that  can  be  read  by  many  CFD  software  packages. 

3.  Theory 

In  this  section  the  background  for  determining  the  homoge¬ 
nised  parameters,  i.e.  porosity,  surface  area,  tortuosity  and  effective 
diffusion  coefficients  will  be  presented. 

3.1.  Porosity 

Once  the  micro-structure  is  reconstructed,  the  determination  of 
an  averaged  porosity  e  is  straightforward: 


where  Vvo  is  the  volume  of  one  pore  voxel,  nvo  the  number  of  pore 
voxels,  and  Ve  the  entire  volume  of  the  considered  domain.  The 
value  of  the  porosity  can  easily  be  determined  by  measurements 
[33]  and  does  not  necessarily  need  the  effort  of  micro-structural 
reconstruction.  However,  it  is  possible  to  retrieve  the  spatial  dis¬ 
tribution  of  the  porosity,  which  is  an  important  information  and  is 
described  in  Section  4.2  and  shown  in  Fig.  3. 

3.2.  Surface  area 

The  reconstruction  algorithm  also  determines  the  specific  sur¬ 
face  area  (surface  area  normalised  by  the  phase  volume)  of  the 
regenerated  mesh.  The  results  would  be  misleading  because  the 
voxel  mesh  overestimates  the  surface  area.  A  rescaling  factor  r  of  11/ 
6  is  proposed  [20],  but  this  only  relates  the  surface  of  sphere  to  its 
circumscribed  cube.  The  reconstructed  mesh  approximates  the 
surface,  as  well  as  the  volume.  Therefore,  a  rescale  factor  for  the 
specific  surface  area  has  to  be  found.  In  the  above  case  the  specific 
surface  area  of  the  sphere  and  the  circumscribed  cube  would  be  the 
same.  Basically,  the  rescale  factor  depends  on  the  voxel  resolution 
and  the  algorithm's  capability  of  recreating  the  structure  with 


Fig.  1.  Left-hand  side  shows  a  single  2-D  image  representing  a  slice  through  the  obtained  3-D  micro-structure  with  X-ray  CT-scans.  The  digitalised  binary  image  is  shown  on  the 
right  hand  side  (white  ->  material,  black  — >  pores).  The  gas-channel  would  be  on  top.  The  lower  parts  show  the  anode  functional  layer  (AFL)  and  the  dense  electrolyte.  The 
evaluation  is  confined  to  the  porous  metal-support. 
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Fig.  2.  a)  Image  stack  obtained  with  X-ray  CT-scans  of  a  sintered  half  cell.  On  the  bottom  the  anode  functional  layer  and  the  gas  tight  electrolyte  can  be  seen.  For  the  reconstruction 
only  the  porous  metal-support  is  evaluated,  b)  Reconstructed  computational  geometry  of  the  pores  with  a  voxel  length  of  1.0  x  10-6  m. 


voxel.  The  relation  between  the  specific  area  of  the  sphere  and  its 
voxelised  cubic  counterpart  is: 


6  6nvo  a 

=  r - “ 

u  snvo 

with 


^  ^vo,A 

The  diameter  of  the  sphere  is  d ,  the  voxel  length  is  s  and  nV0,A 
stands  for  the  number  of  surface  voxel  at  each  side  of  the  voxelised 
sphere.  The  fraction  nvo/nv 0,a  accounts  for  the  algorithm's  capability 
to  reconstruct  the  sphere.  The  mean  uncorrected  specific  surface 
area  for  the  metal-support  is  0.117  pm-1. 
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Fig.  3.  Porosity  distribution  along  the  three  coordinate  axes  for  the  windowing  ge¬ 
ometry  with  a  base  area  of  350  pm  x  350  pm.  The  z-direction  is  the  diffusion  path  from 
the  gas-channel  to  the  anode  functional  layer.  While  the  directions  parallel  to  the  gas- 
channel  show  a  homogeneous  porosity  distribution,  the  diffusion  path  orthogonal  to 
the  gas  channel  has  a  much  lower  porosity  in  the  first  50  to  100  pm  compared  to  the 
rest  of  the  electrode.  The  gas-channel  interface  is  located  at  0  pm  and  the  interface  to 
the  anode  functional  layer  at  350  pm. 


3.3.  Tortuosity 


The  tortuosity  r  is  more  difficult  to  determine.  Many  different 
definitions  of  the  tortuosity  exist  in  literature  [40  .  A  geometrical 
definition  of  the  tortuosity  in  2-D  is  the  ratio  between  the  length  of 
the  actual  or  effective  path  Leff  and  the  straight  and  shortest 
possible  distance  L  between  two  points: 


There  are  attempts  to  determine  the  tortuosity  experimentally 
through  saturation  current  density  measurement  [41],  or  mercury 
intrusion  porosimetry  (MIP)  [42],  although  the  parameter  cannot 
be  measured  directly  [40  .  Computation  of  it,  is  therefore  an 
excellent  alternative.  The  geometrical  tortuosity  can  be  calculated 
by  the  fast  marching  method  43]  or  the  centre  of  mass  calculations 
[17].  Another  possibility  is  to  simulate  a  transport  through  the 
related  phase  (here  pores)  and  analyse  the  impact  of  the  other 
phases  (here  metal  backbone  of  the  support)  on  the  transport 
properties.  In  Ref.  15]  a  tortuosity  was  obtained  by  solving  the 
Laplace  equation  for  the  electrical  potential,  which  provides  the 
ratio  of  effective  conductivity  to  bulk  conductivity  of  the  related 
phase.  In  this  paper  the  tortuosity  is  determined  by  simulating  the 
diffusive  flux  through  the  pores  and  calculating  an  effective  diffu¬ 
sion  coefficient.  In  this  case  the  tortuosity  equals  a  kinematic 
property,  comparable  to  the  relative  average  length  of  the  flow  path 
of  a  fluid  particle  from  one  side  of  a  porous  structure  to  the  other.  Its 
original  introduction  in  SOFC  models  was  also  intended  as  a  kine¬ 
matic  property  13].  The  advantage  is  that  it  is  directly  related  to  the 
effective  transport  properties  of  the  layer  and  thus  can  be  used  for 
performance  simulations. 


3.4.  Effective  diffusion  coefficient 

Fields  law  is  the  simplest  form  to  describe  binary  diffusion  12]. 
In  1-D  the  molar  flux  Na  can  be  calculated  by: 

A  v 

Na  =  -cg,sDab — (5) 

where  Dab  is  the  binary  diffusion  coefficient,  for  example  a  =  H2  and 
b  =  H2O,  Cgas  is  the  gas  concentration,  Axa  is  the  difference  in  the 
mole  fraction,  and  L  is  the  length  between  the  gas-channel  and  the 
anode  functional  layer,  which  is  in  this  case  equivalent  to  the 
thickness  of  the  metal-support.  Consequently,  the  diffusive  fluxes 
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depend  on  their  diffusion  coefficients  and  the  prevailing  driving 
forces  (here:  molar  concentration).  Several  models  are  available  to 
describe  the  binary  diffusion  coefficient  Dat,  see  e.g.  Ref.  [44].  In  the 
presence  of  porous  media,  diffusion  is  impeded  by  the  porosity  and 
the  tortuosity.  This  is  taken  into  account  by  introducing  an  effective 
binary  diffusion  coefficient  D^f  in  porous  media  [44-48  ]: 

jvr°us  =  — CgasDf^ ,  (6) 

where  jv£orous  is  the  superficial  axial  diffusion. 

According  to  Refs.  [33,49]  also  the  underlying  driving  force  is 
influenced  by  the  tortuosity.  The  definition  of  the  average  diffu- 
sional  flux  through  the  pores  is: 

Ax 

Na,p  =  —CgasDafr-z  ,  (7) 

Leff 

where  Leff  is  the  actual  pore  length.  The  superficial  axial  diffusion 
can  be  related  to  the  interstitial  axial  diffusion  Naa  by 

K°r0US  =  e.  (8) 

The  residence  time  t  of  species  a  in  the  porous  medium  can  be 
expressed  either  by  the  interstitial  axial  diffusion: 


L  c 


gas 


N 


a, i 


or  the  diffusional  flux  through  the  pores: 


t  =  L-^.  (10) 

iVa,p 

Combining  equations  (8) — (10)  results  in: 


NaP°rous 

e  L  ' 


(11) 


Inserting  equation  (11)  into  (6)  and  using  the  definition  of  the 
tortuosity  from  equation  (4)  yields  after  rearranging: 


K°mUS  =  -cgas^Dab^L.  (12) 

Comparing  equations  (12)  and  (6)  and  defining  the  scaling  factor 
/as  D^f  =  fDah  it  should  read: 

/= 4-  (13) 

TZ 

as  was  correctly  mentioned  by  Refs.  18,50].  The  effective  diffusion 
coefficient  is  an  important  parameter,  used  to  describe  the  trans¬ 
port  in  homogenised  electrode  models  [51  .  It  is  often  employed  as 
a  fitting  parameter  in  macro-scale  simulations  [13,52-54],  why  the 
misinterpreted  use  of  the  scaling  factor  /  =  e/r  instead  off  =  e/r2 
is  not  crucial.  However,  in  micro-structural  modelling  the  actual 
molar  flux  in  the  reconstructed,  tortuous  geometry  is  computed 
and  compared  with  equation  (6),  which  would  be  the  macroscopic 
approach  to  calculate  the  molar  flux  in  a  homogenised  geometry. 
Thus  it  is  important  to  use  the  correct  formulation  of  the  scaling 
factor  given  in  equation  (13). 

The  procedure  to  calculate  the  scaling  factor  is  presented  in 
Section  4.1.  When  using  the  actual  micro-structure  this  factor  also 
includes  a  correction  for  pockets  and  dead-end  pores.  In  addition 
to  this,  Ref.  [55  explained  that  not  only  the  longer  path  length  has 
to  be  taken  into  account,  but  also  the  existence  of  small  orifices, 
which  impedes  diffusion.  To  sum  up,  the  calculated  scaling  factor 


allows  the  estimation  of  the  tortuosity  of  the  micro-structure,  but 
this  value  does  not  necessarily  equal  the  geometrical  definition  of 
the  tortuosity,  given  in  equation  (4),  for  the  reasons  mentioned 
above. 

4.  Methodology 

4.2.  Calculation  of  the  effective  diffusion  coefficient 

The  solution  of  the  diffusion  transport  equation  provides  the 
distribution  of  the  molar  fraction  x.  This  approach  can  also  be 
applied  for  electrons  or  ions,  as  their  transport  is  described  by  the 
same  formulation.  Considering  steady-state  conditions  and  no 
external  body  forces,  it  can  be  written  for  the  mass  diffusion  as 
follows: 

V(  -  DabVXa)  =  0.  (14) 

Dirichlet  boundary  conditions  x  =  xt,  x  =  xb  are  often  applied  on 
the  top  and  bottom  surfaces  [15,16,18,19,56  ,  while  the  gas-solid 
interface  is  a  Neumann  boundary  (zero-flux).  The  influence  of  the 
boundary  conditions  on  the  top  and  bottom  surfaces  was  investi¬ 
gated  and  Dirichlet  and  Neumann  boundary  conditions  were 
compared  [20  .  It  was  shown,  that  different  results  were  obtained 
and  an  intermediate  boundary  condition  (fixed  concentration  on 
top  and  fixed  flux  on  bottom)  was  proposed  as  best  practice.  In  all 
these  works  the  remaining  boundaries  of  the  computational 
domain  are  typically  set  to  a  zero-flux  Neumann  condition  Vx  =  0, 
as  this  is  often  applied  to  represent  symmetry.  It  was  discussed  in 
the  Introduction,  that  this  models  reality  only  in  an  approximate 
way,  as  this  effectively  introduces  a  lot  of  dead-end  pores,  or  in 
other  words  blocks  tortuous  paths  involving  diffusion  in  and  out  of 
the  domain  in  the  vicinity  of  the  lateral  boundary. 


Fig.  4.  Pore  volume  of  the  computational  geometry  with  lateral  boundaries  shown  in 
dark  blue  ->  zero  gradient  and  red  ->■  periodic  boundary.  Periodic  boundary  condi¬ 
tions  are  applied  wherever  possible.  The  size  of  the  geometry  is 
300  pm  x  300  pm  x  350  pm  and  the  diffusion  path  is  from  top  to  bottom.  (For 
interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to 
the  web  version  of  this  article.) 


G.  Reiss  et  al.  /  Journal  of  Power  Sources  273  (2015)  1006-1015 


1011 


A  more  realistic  approach  is  to  use  so-called  periodic  boundaries 
for  this  lateral  boundary  to  the  highest  feasible  extent,  as  is  done  in 
many  different  fields  [36,57-59].  This  means  that  the  mass,  which 
leaves  the  computational  geometry  at  one  side,  enters  at  the 
opposite.  This  would  also  be  the  case  in  reality,  where  the  pores  will 
leave  and  enter  the  examined  geometry.  In  Fig.  4  the  computational 
geometry  is  shown.  The  diffusion  path  is  from  top  to  bottom  and 
the  lateral  boundary  are  shown.  For  boundary  regions,  where  a 
porous  path  is  continued  on  the  opposite  side  of  the  considered 
domain,  a  periodic  boundary  condition  is  set.  These  boundary  re¬ 
gions  are  coloured  red.  The  blue  coloured  boundary  regions  sig¬ 
nifies  porous  paths  with  no  counterpart  on  the  opposite  side,  why 
they  are  effectively  blocked  in  this  simulation.  These  are  coloured 
blue.  This  new  approach  will  lead  to  a  better  prediction  of  the 
scaling  factor  because  the  zero-flux  Neumann  condition  can  lead  to 
an  overestimation  of  the  scaling  factor,  since  many  of  the  prevailing 
pores  are  then  supposed  to  be  dead-end  pores. 

The  diffusion  equation  (14)  is  numerically  solved  by  integrating 
over  the  computational  cell  and  applying  the  Gauss-theorem.  The 
surface  integral  is  then  discretized  using  the  finite  volume  method 
(FVM)  based  on  the  open-source  CFD  toolkit  OpenFOAM  [60]: 

J  V(DabVXa)dV  =  JfSf(DabVXa)  =  0,  (15) 


(18) 


4.2.  Calculation  of  a  3-D  porosity  distribution 

In  Section  3.1  the  calculation  of  an  averaged  porosity  is 
described,  which  can  be  used  in  cell  and  stack  simulation.  In  order 
to  justify  the  efforts  of  a  reconstruction  a  MATLAB  algorithm  was 
programmed,  which  analyses  the  micro-structure  and  computes 
i.e.  the  porosity  distribution,  the  mean  pore  diameter,  the  pore 
number,  the  pore  area  size  and  orientation  along  all  directions.  This 
gives  a  realistic  picture  of  the  inhomogeneity  of  the  micro-structure 
leading  to  anisotropic  transport  properties.  The  basic  procedure  of 
the  algorithm  is  as  follows: 

1  The  coordinates  of  each  computational  cell  serve  as  an  input 
parameter  for  the  algorithm. 

2  One  coordinate  axis  can  be  chosen.  The  sample  is  subdivided 
into  slices  orthogonal  to  this  axis  (see  left  hand  side  Fig.  5). 

3  Every  slice  is  processed  using  the  standard  MATLAB  Image  Pro¬ 
cessing  Toolbox,  which  distinguishes  spatially  separated  pore 
cross  sections  (see  right  hand  side  Fig.  5). 

4  These  cross-sections  are  analysed,  and  the  results  can  be  plotted 
along  the  previous  defined  axis. 


where  Vp  is  the  control  volume  of  one  computational  cell,  the 
subscript /denotes  a  face  of  this  control  volume,  and  Sf  is  the  sur¬ 
face  normal  vector  of  the  face.  Once  the  distribution  of  x  is  calcu¬ 
lated  the  averaged  local  flux  f  across  the  boundary  surface  can  be 
calculated: 


Exemplary  results  of  the  porosity  distribution  for  different 
sample  positions  are  plotted  in  Fig.  3.  The  distribution  of  the  mean 
pore  diameter  is  shown  in  Fig.  6,  and  the  pore  number  is  presented 

in  Fig.  10. 


-DabVxandS, 


(16) 


and  hence  the  effective  diffusion  coefficient  can  be  determined: 


s 


-DabVxandS 


(¥)* 


(17) 


The  denominator  in  equation  (17)  represents  the  molar  flux 
calculated  with  equation  (6)  in  the  homogenised  geometry  across 
the  area  A.  The  difference  in  the  mole  fraction  A*a  is  calculated  as 
the  difference  of  the  weighted  averaged  mole  fraction  at  the  top 
and  bottom  surfaces.  Here  the  same  approach  as  Ref.  20  proposed, 
is  pursued  and  different  boundary  conditions  (Dirichlet,  Neumann) 
are  applied  to  the  top  and  bottom  boundaries.  A  similar  tendency  is 
observed: 


1  Two  fixed  mole  fractions  on  the  top  and  bottom  boundaries 
(Dirichlet)  result  in  the  lowest  tortuosity  value. 

2  Two  fixed  fluxes  on  the  top  and  bottom  boundaries  (Neumann) 
yield  the  highest  tortuosity  value. 

3  The  tortuosity  value  obtained  with  a  fixed  mole  fraction  on  the 
top  and  a  defined  flux  on  the  bottom  boundary  lies  between  the 
above  approaches. 

In  this  paper  the  last  approach  is  used  as  this  resembles  the 
nature  of  an  SOFC  most  adequately,  see  also  Ref.  [20  . 

With  the  effective  diffusion  coefficient  and  the  definition  of  the 
scaling  factor/in  equation  (13)  the  tortuosity  of  the  micro-structure 
can  be  determined: 


4.3.  Representative  reconstructed  geometries 

The  procedure  for  calculating  averaged  parameters  out  of  a 
reconstructed  micro-structure  is  limited  by  computational  capacity. 
The  larger  the  investigated  volume,  the  more  computational  cells 
are  needed  to  mesh  the  domain,  the  longer  are  the  calculation  times, 
but  it  is  more  representative  for  the  entire  cell.  When  dealing  with 
an  inhomogeneous  material  the  results  of  the  averaging  processes 
depend  on  the  position,  where  the  sample  is  taken  and  on  the 
sample  size  [33  .  In  order  to  decide  whether  or  not  the  sample  size  is 
large  enough  a  stochastic  equivalent  representative  volume  element 
(SERVE)  is  determined  [15,16,20],  by  the  so-called  windowing 
technique.  In  windowing  the  volume  of  the  investigated  geometry  is 
increased  until  e.g.  the  change  of  the  desired  homogenised  param¬ 
eter  (e.g.  porosity)  drops  below  a  desired  threshold.  However,  the 
analysis  is  limited  to  one  location  of  this  sample,  which  was  taken 
out  of  one  place  of  a  larger  cell.  This  approach  is  applicable,  if  it  can 
be  assumed  that  the  micro-structure  is  more  or  less  homogeneous. 
In  order  to  get  a  more  representative  result,  it  is  important  to  use 
reconstructed  geometries  from  different  positions  of  the  sample,  as 
done  in  this  work.  The  importance  of  this  can  be  understood  when 
looking  at  Fig.  3.  The  inhomogeneous  nature  of  the  metal-support 
can  be  seen  in  the  huge  difference  of  the  porosity  in  the  first 
100  pm  compared  to  the  remaining  part.  Determining  a  SERVE  and 
starting  from  the  wrong  position  may  lead  to  wrong  results. 
Therefore,  the  entire  height  of  350  pm  is  considered  in  this  paper 
and  the  base  area  of  the  reconstructed  geometries  is  increased, 
while  the  centre  point  stays  the  same  for  each  geometry. 

5.  Results 

The  transport  parameters  are  evaluated  using  a  3-D  recon¬ 
struction  of  a  metal-support  with  subsequent  numerical 
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Fig.  5.  Computational  geometry  is  cut  into  slices  (left  hand  side).  These  slices  are  analysed  with  MATLAB  and  their  size,  orientation,  shape,  area  and  diameter  is  determined.  Right 
hand  side  shows  an  evaluated  sample  slice,  where  connected  pores  are  highlighted  in  different  colours.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the 
reader  is  referred  to  the  web  version  of  this  article.) 
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Fig.  6.  Mean  pore  diameter  distribution  along  the  diffusion  path  for  two  different 
positions  in  the  sample.  The  gas-channel  interface  will  be  at  0  pm  and  the  interface  to 
the  anode  functional  layer  at  350  pm. 


computation.  The  obtained  homogenised  parameters  for  the 
porous  metal-support  are  shown  in  Table  1  and  the  different 
sample  positions  are  depicted  in  Fig.  7.  The  homogenised  mean 
porosity  for  five  different  sample  positions  is  0.200  and  the  mean 
equivalent  diameter  is  18.1  pm.  The  value  for  the  homogenised 
computed  tortuosity  is  4.3  calculated  with  equation  (18).  It  was 
discussed  in  Section  3.4  and  stressed  again  here,  that  the  calculated 
tortuosity  from  the  micro-structural  computation  also  includes 
effects  like  dead-end  pores,  orifices  and  pockets  and  does  not 
necessarily  equal  the  geometrical  definition  of  the  tortuosity  given 
in  equation  (4).  It  has  to  be  noted,  that  these  values  correspond  to 
an  uncorroded  micro-structure  and  may  decrease  during  operation 
and  expected  corrosion  of  the  metal-support.  This  issue  is  subject 
to  further  development  of  a  corrosion  model. 

The  computed  scaling  factors  are  compared  to  measurement 
results  already  published  for  a  similar  cell  design  [2  .  There  are  only 
minor  differences  in  porosity  and  pore-size  distribution  of  the 
metal-support  between  the  simulated  and  the  tested  cell.  Based  on 
electro-chemical  impedance  spectroscopy  and  a  subsequent 
impedance  data  analysis  by  the  distribution  of  relaxation  times,  the 
individual  loss  mechanisms  could  be  identified.  The  area  specific 
resistance,  which  is  related  to  gas  diffusion  processes  in  the  sub¬ 
strate,  is  extracted  and  fitted  with  a  model  proposed  by  Ref.  14]. 
The  fitting  result  for  the  scaling  factor  /,  which  describes  the 
diffusion  process  in  the  metal-support,  is  1.32  x  10-2.  The  com¬ 
parison  of  the  results  is  summarised  in  fable  2,  where  the  statistical 


Table  1 

Comparison  of  micro-structural  parameters,  obtained  from  3-D  reconstruction/ 
modelling,  for  five  different  computational  geometries,  with  a  base  area  size  of 
200  pm  x  200  pm.  The  subscript  horn  refers  to  the  overall  homogenised  parameter 
calculated  with  a  height  of  350  pm.  The  dense-parameters  are  evaluated  only  in  the 
first  100  pm,  where  a  significant  lower  porosity  was  found  (see  Fig.  3).  Parameters 
with  the  subscript  porous  are  computed  in  the  remaining  region  (from  100  to 
350  pm).  The  principal  positions  of  the  geometries  are  shown  in  Fig.  7,  and  Fig.  8 
depicts  the  positions  of  the  different  layers. 


Mean  value 

Standard  deviation 

^hom 

0.200 

0.009 

Scaling  factor /hom 

1.19  x  102 

4.8  x  10  3 

Thom 

4.3 

0.81 

Pore  diameter  (pm) 

18.1 

0.51 

^dense 

0.147 

0.006 

Scaling  factor  /dense 

0.84  x  10  2 

4.03  x  10  3 

^porous 

0.219 

0.013 

Scaling  factor  /porous 

1.94  x  10  2 

5.7  x  10~3 
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Another  interesting  finding,  which  can  be  seen  in  ^able  2  is  that 
the  difference  between  the  zero  gradient  and  periodic  boundary 
conditions  is  about  46.9%  and  that  the  more  accurate  periodic 
approach  shows  a  better  agreement  with  the  measurement  data. 
This  supports  and  proves  the  assumption  that  the  incorrect  appli¬ 
cation  of  zero  gradient  boundary  conditions  at  the  lateral  walls 
introduces  additional  dead-end  pores  which  reduce  the  entire  mass 
transport  through  the  micro-structure  and  leads  to  a  smaller  the 
scaling  factor. 

One  explanation  for  the  low  scaling  factor  is,  that  this  could  be 
due  to  the  dense  first  100  pm  from  the  gas-channel  to  the  anode 
functional  layer,  where  also  the  number  of  pores  is  much  lower 
compared  to  the  rest  of  the  metal-support,  as  shown  in  Fig.  10. 
Therefore,  a  more  thorough  analysis  of  the  metal-support  micro¬ 
structure  is  conducted.  In  addition  to  the  homogenised  transport 
parameters  for  the  whole  micro-structure,  the  denser  first  100  pm, 
as  well  as,  the  remaining  more  porous  250  pm  are  assessed  sepa¬ 
rately,  see  Fig.  8  for  a  sketch  of  the  layers.  A  comparison  is  shown  in 
Table  1,  again  with  the  statistical  values  of  five  different  computa¬ 
tional  geometries.  Going  from  the  gas-channel  to  the  anode  func¬ 
tional  layer  the  first  100  pm  of  the  metal-support  have  an  averaged 
porosity  of  0.147  and  a  scaling  factor  of  0.84  x  1CT2.  The  region  from 
100  to  350  pm  has  an  averaged  porosity  of  0.219  and  a  scaling  factor 
of  1.94  x  1CT2.  The  difference  in  the  scaling  factors  is  about  130%, 
which  allows  to  draw  the  conclusion,  that  the  results  of  a  cell 
model,  calculated  with  two  different  layers  will  be  more  accurate 
than  those  using  an  overall  homogenised  approach. 

Another  problem  addressed,  is  the  question  of  how  to  find  a 
stochastic  equivalent  representative  volume  element  (SERVE).  The 
state-of-the-art  approach  is  explained  and  its  drawbacks  can  be 
assessed  by  this  inhomogeneous  micro-structure.  As  shown  in 
Fig.  3  the  porosity  in  the  first  100  pm,  from  the  gas-channel  to  the 
anode  functional  layer,  is  much  lower  compared  to  the  rest  of  the 


Fig.  7.  Position  of  the  arbitrarily  chosen  samples.  It  is  ensured  that  the  sample  posi¬ 
tions  don't  overlap  and  a  certain  distance  is  kept  to  the  edges. 

values  of  five  different,  reconstructed  geometries  are  presented. 
The  positions  of  these  geometries  are  shown  in  Fig.  7.  Comparing 
the  computed  and  measured/fitted  scaling  factor  it  can  be  seen  that 
the  values  are  in  the  same  order  of  magnitude  and  that  the  sug¬ 
gested  approach  with  periodic  boundary  conditions  yields  a  value 
of  1.19  x  10-2,  which  deviates  only  by  11%  from  the  measurements. 
Given  the  differences  of  the  two  studied  cells  this  is  an  expected 
deviation.  Furthermore,  the  scaling  factor  in  Table  2  of  the  metal- 
support  is  much  smaller,  by  a  factor  of  20,  compared  to  values, 
which  are  commonly  employed  in  macro-scale  simulations  and  are 
in  the  range  of  10-25  x  1CT2  [47,48,50,54,61,62].  However,  it  has  to 
be  mentioned,  that  in  these  works  the  scaling  factor  was  used  as  a 
fitting  parameter  for  the  applied  simulation  models,  and  was  not 
determined  either  by  measurement/fitting  or  3-D  reconstruction/ 
modelling. 


Table  2 

Comparison  of  the  mean  values  and  standard  deviations  of  scaling  factors  for 
different  boundary  conditions  calculated  on  the  same  computational  geometries  to 
the  value  obtained  from  impedance  measurements  which  is  1.32  x  10~2.  The  base 
area  is  in  all  cases  200  pm  x  200  pm  and  the  height  is  350  pm. 


Boundary  condition  (BC) 

Scaling  factor  (xlO-2) 

Mean  value  Standard  deviation 

Diff  (%) 

Measurement 

1.32 

Zero  gradient 

0.81 

0.41 

63.0 

Periodic 

1.19 

0.477 

10.9 

gas-channel 


Fig.  8.  Principal  depiction  of  the  metal-support  with  the  100  pm  thick  dense  upper 
layer  (referred  to  as  dense)  and  the  more  porous  250  pm  of  the  remaining  metal- 
support  (referred  to  as  porous).  Homogenised  parameters  (referred  to  as  horn)  are 
calculated  for  the  entire  height  of  350  pm. 
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Table  3 

Porosity  values,  scaling  factor  and  number  of  computational  cells  for  a  geometry 
height  of  350  pm  and  increasing  base  area.  The  geometry  with  a  base  area  of 
50  pm  x  50  pm  has  not  got  a  continuous  pore  from  the  gas-channel  to  the  anode 
functional  layer  interface. 


Area  Odense 

Porosity  (-) 

Scaling 

factor  (-)  (xl0~2) 

Number  of 
computational  cells 

50  x  50 

0.193 

— 

0.2  x  106 

100  x  100 

0.196 

2.11 

0.7  x  106 

150  x  150 

0.206 

1.90 

1.6  x  106 

200  x  200 

0.206 

1.99 

2.9  x  106 

250  x  250 

0.192 

1.69 

4.2  x  106 

300  x  300 

0.193 

1.35 

6.1  x  106 

350  x  350 

0.209 

1.73 

9.0  x  106 

micro-structure.  This  may  be  due  to  the  sintering  process,  and  is 
subject  to  further  investigation.  When  a  SERVE  is  determined  the 
volume  of  the  reconstructed  geometry  is  increased  till  the  param¬ 
eters  remain  constant  within  a  certain  tolerance.  Nevertheless, 
increasing  the  SERVE  and  starting  from  the  anode  functional  layer,  a 
SERVE  might  be  found  in  the  metal-support  micro-structure  before 
the  dense  region  near  the  gas-channel  is  reached.  This  means  that 
without  the  knowledge  about  the  (in-)homogeneity  of  the  micro¬ 
structure  the  starting  point  of  the  SERVE  investigation  may  lead 
to  a  wrong  assessment  of  the  micro-structure.  Therefore,  it  was 
attempted  to  use  the  entire  height  of  the  metal-support  for  the 
determination  of  a  SERVE  and  the  assessment  of  transport  pa¬ 
rameters  in  this  paper.  Starting  with  a  base  area  of  50  pm  x  50  pm, 
the  area  is  gradually  increased  and  the  micro-structural  parameters 
determined.  The  porosity  value  and  the  scaling  factor  along  with 
the  number  of  computational  cells  are  listed  in  Table  3.  Although 
the  volume  of  the  reconstructed  geometry  is  increased,  the  porosity 
and  the  tortuosity  vary  and  no  clear  convergence  of  the  parameters 
can  be  found.  In  Fig.  9  the  porosity  distribution  is  plotted  along  the 
diffusion  path  for  increased  base  areas.  Here  again  the  deviation  of 
the  porosity  decreases,  but  the  volume  is  still  not  sufficiently  large, 
although  it  exceeds  the  size  of  SERVES  reported  in  literature  by 
orders  of  magnitude  [16-18,20  .  Furthermore,  the  rather  high  de¬ 
viations  of  the  scaling  factors  in  Table  2  indicate  that  the  metal- 
support  is  inhomogeneous  and  the  results  of  the  computational 
evaluation  of  the  micro-structure  highly  depends  on  the  chosen 


Fig.  9.  Porosity  distribution  for  reconstructed  geometries  with  the  same  centre  point 
and  increasing  base  areas.  Although  the  base  area  is  increased  the  porosity  distribution 
still  changes  and  no  SERVE  is  obtained.  The  gas-channel  interface  will  be  at  0  pm  and 
the  interface  to  the  anode  functional  layer  at  350  pm. 
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Fig.  10.  Number  of  pores  per  m2  along  the  diffusion  path.  In  the  first  100  pm  there  are 
less  pores  compared  to  the  rest  of  the  metal-support.  The  gas-channel  interface  will  be 
at  0  pm  and  the  interface  to  the  anode  functional  layer  at  350  pm.  The  base  area  is 
200  pm  x  200  pm. 


location.  This  confirms  the  approach  of  this  paper  to  use  several 
different  sample  locations  and  average  the  retrieved  micro- 
structural  parameters,  as  well  as,  to  distinguish  different  layers 
and  to  determine  the  individual  transport  parameters  separately. 


6.  Conclusion 

In  this  paper  the  computational  framework  for  the  investigation 
of  homogenised  parameters  of  a  porous  metal-support  for  metal- 
supported  SOFC  has  been  presented.  A  more  realistic  set  of 
boundary  conditions  for  homogenisation  was  introduced.  It  was 
shown  that  for  the  considered  micro-structure  a  significant  error  is 
committed  by  applying  the  commonly  used  zero  flux  boundary 
condition  on  the  lateral  boundaries,  and  that  the  usage  of  the  more 
accurate  periodic  boundary  condition  to  the  highest  feasible  extent 
is  recommended.  The  values  of  the  homogenised  scaling  factor  for 
the  effective  diffusion  coefficient  for  the  metallic  support  was 
determined  to  be  approximately  1.2  x  1CT2.  The  results  are 
compared  to  impedance  spectroscopy  measurement  and  showed 
good  agreement.  In  order  to  assess  the  inhomogeneity  of  the  micro¬ 
structure  the  porosity  distribution,  as  well  as,  the  transport  pa¬ 
rameters  of  two  identified,  different  dense  layers  of  the  metal- 
support  were  evaluated.  The  presented  model  is  a  first  step  to¬ 
ward  a  study  of  the  influence  of  corrosion  on  the  effective  transport 
parameters. 
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Nomenclature 

e  porosity  (-) 

r  tortuosity  (-) 

r2  tortuosity  factor  (-) 

0  flux  e.g.  molar  flux  (kmol  rrT2  s'1) 

Na  molar  flux  of  species  a  (kmol  m'2  s'1) 

j^poious  mo|ar  fjux  0f  Species  a  in  porous  medium(kmol  m'2  s'1) 

Na,i  axial  interstitial  molar  flux  of  species  a  in  porous  medium 
(kmol  m~2  s'1) 

Na  p  average  diffusional  molar  flux  of  species  a  through  the 
pores  of  a  porous  medium  (kmol  m'2  s'1) 

Sf  surface  normal  vector  of  face  (m2) 

A  base  area  of  the  homogenised  geometry  (m2) 

cgas  concentration  (kmol  m  3) 

d  diameter  of  a  sphere  (m) 

Dab  binary  diffusion  coefficient  (m2  s'1) 

D^f  effective  binary  diffusion  coefficient  (m2  s'1) 

f  scaling  factor  (-) 

L  length  (m) 

Leff  effective  length  (m) 

nvo, a  number  of  voxel  at  the  surface  (-) 

nvo  number  of  voxel  (-) 

r  rescaling  factor  for  specific  surface  area  (-) 

s  length  of  one  voxel  (m) 

t  residence  time  (s) 

Ve  volume  of  computational  domain  (m3) 

Up  control  volume  of  one  computational  cell  (m3) 

Uvo  volume  of  one  voxel  (m3) 

xa  mole  fraction  of  species  a  (-) 

References 

[1]  P.  Blennow,  J.  Hjelm,  T.  Klemenso,  A.H.  Persson,  S.  Ramousse,  M.  Mogensen, 
Fuel  Cells  11  (2011)  661-668. 

[2]  P.  Blennow,  J.  Hjelm,  T.  Klemenso,  S.  Ramousse,  A.  Kromp,  A.  Leonide, 
A.  Weber,  J.  Power  Sources  196  (2011)  7117—7125,  9th  European  SOFC  Forum, 
Lucerne,  Switzerland,  2010. 

[3]  M.C.  Tucker,  G.Y.  Lau,  C.P.  Jacobson,  L.C.  Dejonghe,  S.J.  Visco,  J.  Power  Sources 
171  (2007)  477-482. 

[4]  M.C.  Tucker,  G.Y.  Lau,  C.P.  Jacobson,  L.C.  Dejonghe,  S.J.  Visco,  J.  Power  Sources 
175  (2008)  447-451. 

[5]  M.C.  Tucker,  J.  Power  Sources  195  (2010)  4570-4582. 

[6]  Y.B.  Matus,  L.C.D.  Jonghe,  C.P.  Jacobson,  S.J.  Visco,  Solid  State  Ionics  176  (2005) 
443-449. 

[7]  T.  Franco,  M.  Haydn,  R.  Miicke,  A.  Weber,  M.  Riittinger,  O.  Biichler, 
S.  Uhlenbruck,  N.H.  Menzler,  A.  Venskutonis,  L.S.  Sigl,  ECS  Trans.  35  (2011) 
343-349. 

[8]  T.  Franco,  R.  Miicke,  A.  Weber,  M.  Haydn,  M.  Riittinger,  N.H.  Menzler,  A. 
Venskutonis,  L.S.  Sigl,  H.P.  Buchkremer,  in:  F.  Lefebbvre-Joud,  Lefe  (Eds.), 
Proceedings  of  the  10th  European  Solid  Oxide  Fuel  Cell  Forum,  Lucerne, 
Switzerland,  A0906  49— A0906  61. 

[9]  T.  Franco,  M.  Haydn,  A.  Weber,  W.  Schafbauer,  L.  Blum,  U.  Packbier, 
D.  Roehrens,  N.H.  Menzler,  J.  Rechberger,  A.  Venskutonis,  L.S.  Sigl, 
H.P.  Buchkremer,  ECS  Trans.  57  (2013)  471-480. 

[10]  J.  Nielsen,  T.  Klemenso,  P.  Blennow,  J.  Power  Sources  219  (2012)  305—316. 

[11]  B.  Haberman,  J.  Young,  Int.  J.  Heat  Mass  Transf.  47  (2004)  3617—3629. 

[12]  R.  Bove,  S.  Ubertini,  Modeling  Solid  Oxide  Fuel  Cells:  Methods,  Procedures  and 
Techniques,  Springer,  2008. 

[13]  M.  Andersson,  J.  Yuan,  B.  Sunden,  Appl.  Energy  87  (2010)  1461—1476. 

[14]  A.  Leonide,  V.  Sonn,  A.  Weber,  E.  Ivers-Tiffee,  J.  Electrochem.  Soc.  155  (2008) 
B36-B41. 

[15]  J.  Joos,  T.  Carraro,  A.  Weber,  E.  Ivers-Tiffe,  J.  Power  Sources  196  (2011) 
7302-7307. 

[16]  J.  Joos,  M.  Ender,  T.  Carraro,  A.  Weber,  E.  Ivers-Tiffe,  Electrochim.  Acta  82 
(2012)  268-276. 


[17]  D.  Gostovic,  J.R.  Smith,  D.P.  Kundinger,  K.S.  Jones,  E.D.  Wachsman,  Electro¬ 
chem.  Solid-State  Lett.  10  (2007)  B214-B217. 

[18]  N.  Vivet,  S.  Chupin,  E.  Estrade,  T.  Piquero,  P.  Pommier,  D.  Rochais,  E.  Bruneton, 
J.  Power  Sources  196  (2011)  7541-7549. 

[19]  N.S.K.  Gunda,  S.K.  Mitra,  J.  Phys.  Conf.  Ser.  362  (2012). 

[20]  J.  Laurencin,  R.  Quey,  G.  Delette,  H.  Suhonen,  P.  Cloetens,  P.  Bleuet,  J.  Power 
Sources  198  (2012)  182-189. 

[21]  J.R.  Izzo,  A.S.  Joshi,  K.N.  Grew,  W.K.S.  Chiu,  A.  Tkachuk,  S.H.  Wang,  W.  Yun, 
J.  Electrochem.  Soc.  155  (2008)  B504. 

[22]  P.  Jorgensen,  K.  Hansen,  R.  Larsen,  J.  Bowen,  J.  Power  Sources  195  (2010) 
8168-8176. 

[23]  G.J.  Nelson,  W.M.  Harris,  J.J.  Lombardo,  J.R.  Izzo  Jr.,  W.K.  Chiu,  P.  Tanasini, 
M.  Cantoni,  J.  Van  herle,  C.  Comninellis,  J.C.  Andrews,  Y.  Liu,  P.  Pianetta, 
Y.S.  Chu,  Electrochem.  Commun.  13  (2011)  586—589. 

[24]  P.  Shearing,  J.  Gelb,  N.  Brandon,  J.  Eur.  Ceram.  Soc.  30  (2010)  1809—1814. 

[25]  G.  Borgefors,  I.  Nystrm,  Pattern  Recognit.  Lett.  18  (1997)  465—471. 

[26]  J.  Wang,  Y.  Tan,  Pattern  Recognit.  46  (2013)  230—242. 

[27]  W.-P.  Choi,  K.-M.  Lam,  W.-C.  Siu,  Pattern  Recognit.  36  (2003)  721-729. 

[28]  C.L.Y.  Yeong,  S.  Torquato,  Phys.  Rev.  E  57  (1998)  495-506. 

[29]  S.  Molin,  B.  Kusz,  M.  Gazda,  P.  Jasinski,  J.  Power  Sources  181  (2008)  3137. 

[30]  S.  Molin,  M.  Gazda,  B.  Kusz,  P.  Jasinski,  J.  Eur.  Ceram.  Soc.  29  (2009)  757762. 

[31]  S.  Molin,  M.  Gazda,  P.  Jasinski,  Solid  State  Ionics  181  (2010)  12141220. 

[32]  S.  Molin,  M.  Gazda,  P.  Jasinski,  J.  Eur.  Ceram.  Soc.  31  (2011)  27072710. 

[33]  J.  Bear,  Dynamics  of  Fluids  in  Porous  Media,  1988.  Dover,  New  York. 

[34]  H.  Yanagihara,  W.  Brandstatter,  N.  Ohashi,  B.  Gschaider,  J.  Leixnering, 
I.  Stankovic,  Top.  Catal.  52  (2009)  1842—1846. 

[35]  E.M.  Borla,  G.  Nicol,  W.  Brandstatter,  M.  Klug,  New  Concepts  and  Technologies 
for  Integrated  Diesel  Exhaust  Gas  Aftertreatment  Systems,  SAE  International, 
Warrendale,  PA,  2011.  SAE  Technical  Paper  2011-24-0184. 

[36]  A.  Matine,  N.  Boyard,  P.  Cartraud,  G.  Legrain,  Y.  Jarny,  Int.  J.  Heat  Mass  Transf. 
62  (2013)  586-603. 

[37]  T.  Klemenso,  J.  Nielsen,  P.  Blennow,  S.H.  Persson,  T.  Stegk,  B.H.  Christensen, 
S.  Snderby,  J.  Power  Sources  196  (2011)  9459—9466. 

[38]  A.  Kromp,  J.  Nielsen,  P.  Blennow,  T.  Klemenso,  A.  Weber,  Fuel  Cells  13  (2013) 
598604. 

[39]  J.  Leixnering,  B.  Gschaider,  W.  Brandsttter,  R.  Bouwman,  in:  OpenFOAM  In¬ 
ternational  Conference,  Old  Windsor,  UK. 

[40]  L.  Shen,  Z.  Chen,  Chem.  Eng.  Sci.  62  (2007)  3748-3755. 

[41]  C.-L.  Tsai,  V.H.  Schmidt,  J.  Power  Sources  196  (2011)  692—699. 

[42]  K.N.  Grew,  A. A.  Peracchio,  A.S.  Joshi,  J.R.  Izzo  Jr.,  W.K.  Chiu,  J.  Power  Sources 
195  (2010)  7930-7942. 

[43]  P.S.  Jorgensen,  K.V.  Hansen,  R.  Larsen,  J.R.  Bowen,  J.  Microsc.  244  (2011) 
45-58.  PMID:  21631501. 

[44]  V.  Novaresio,  M.  Garcia-Camprubi,  S.  Izquierdo,  P.  Asinari,  N.  Fueyo,  Comput. 
Phys.  Commun.  183  (2012)  125—146. 

[45]  M.  Garcia-Camprubi,  A.  Sanchez-Insa,  N.  Fueyo,  Chem.  Eng.  Sci.  65  (2010) 
1668-1677. 

[46]  L.  Pisani,  Int.  J.  Heat  Mass  Transf.  51  (2008)  650-660. 

[47]  W.G.  Bessler,  S.  Gewies,  M.  Vogler,  Electrochim.  Acta  53  (2007)  1782—1800. 

[48]  T.X.  Ho,  P.  Kosinski,  A.C.  Hoffmann,  A.  Vik,  Chem.  Eng.  Sci.  63  (2008) 
5356-5365. 

[49]  N.  Epstein,  Chem.  Eng.  Sci.  44  (1989)  777-779. 

[50]  S.  Wang,  W.M.  Worek,  W.J.  Minkowycz,  Int.  J.  Heat  Mass  Transf.  55  (2012) 
3933-3945. 

[51]  S.B.  Adler,  J.A.  Lane,  B.C.H.  Steele,  J.  Electrochem.  Soc.  143  (1996)  3554—3564. 

[52]  F.N.  Cayan,  S.R.  Pakalapati,  F.  Elizalde-Blancas,  I.  Celik,  J.  Power  Sources  192 
(2009)  467-474. 

[53]  Y.  Vural,  L.  Ma,  D.B.  Ingham,  M.  Pourkashanian,  J.  Power  Sources  195  (2010) 
4893-4904. 

[54]  R.  Suwanwarangkul,  E.  Croiset,  M.W.  Fowler,  P.L.  Douglas,  E.  Entchev, 
M.A.  Douglas,  J.  Power  Sources  122  (2003)  9—18. 

[55]  T.G.  Sherwood,  R.L.  Pigford,  C.R.  Wilke,  in:  B.  Clark,  J.  Maisel  (Eds.),  Mass 
Transfer,  McGraw-Hill  Inc.,  1975,  pp.  39—43. 

[56]  H.  Choi,  A.  Berson,  J.G.  Pharoah,  S.B.  Beale,  Part  A  J.  Power  Energy  225  (2011) 
183-197. 

[57]  S.  Jaganathan,  B.  Pourdeyhimi,  H.V.  Tafreshi,  J.  Eng.  Fibers  Fabr.  3  (2008). 

[58]  A.  Nabovati,  E.W.  Llewellin,  A.C.M.  Sousa,  Compos.  Part  A  Appl.  Sci.  Manuf.  40 
(2009)  860-869. 

[59]  V.D.  Nguyen,  E.  Bchet,  C.  Geuzaine,  L.  Noels,  Comput.  Mater.  Sci.  55  (2011) 
390-406. 

[60]  H.G.  Weller,  G.  Tabor,  H.  Jasak,  C.  Fureby,  Comput.  Phys.  12  (1998)  620—631. 

[61]  S.H.  Chan,  K.A.  Khor,  Z.T.  Xia,  J.  Power  Sources  93  (2001)  130-140. 

[62]  K.  Tseronis,  I.K.  Kookos,  C.  Theodoropoulos,  Chem.  Eng.  Sci.  63  (2008) 
5626-5638. 


